library(tidyverse)
library(haven)
library(arm)

setwd("~/Dropbox/GuessNyhanReifler/DART0023/PNAS/post_accept/replication files/US data/")
dat <- read_dta("US_data_clean_w1_headline.DTA")


mod1 <- lmer(accuracy ~ tips + (tips | dv), data = dat) # random intercept + random slope - all headlines

# all headlines

ggdat1 <- ranef(mod1)$dv %>% dplyr::select(1) %>% rownames_to_column() %>% rename(estimate = "(Intercept)",
                                                                                  headline = rowname) %>% mutate(type = "intercept")
ggdat2 <- ranef(mod1)$dv %>% dplyr::select(2) %>% rownames_to_column() %>% rename(estimate = tips,
                                                                                  headline = rowname) %>% mutate(type = "slope")
ggdat <- bind_rows(ggdat1, ggdat2)
ggdat$SE <- c(se.ranef(mod1)$dv[,1], se.ranef(mod1)$dv[,2])

ggdat$pro <- rep(rep(c(rep("D", 4), rep("R", 4)), 2), 2)
ggdat$source <- rep(c(rep(c(rep("Hyper", 2), rep("False", 2)), 2), rep(c(rep("Mainstream (low prominence)", 2), rep("Mainstream", 2)), 2)), 2)
ggdat$headline <- paste0(ggdat$source, " pro-", ggdat$pro)
ggdat$headline <- paste0(ggdat$headline, " ", rep(1:2, nrow(ggdat)/2))

library(cowplot)

g <- ggplot(filter(ggdat, type == "intercept"), aes(headline, estimate)) + 
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) + 
  geom_pointrange(aes(y = estimate, 
                      ymin = estimate - SE*1.96,
                      ymax = estimate + SE*1.96)) +
  xlab("") + ylab("") + ggtitle("Baseline perceived accuracy") + ylim(-.7, .7) +
  coord_flip()

gg_intercept <- g + theme_minimal()

g <- ggplot(filter(ggdat, type == "slope"), aes(headline, estimate)) + 
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) + 
  geom_pointrange(aes(y = estimate, 
                      ymin = estimate - SE*1.96,
                      ymax = estimate + SE*1.96)) +
  xlab("") + ylab("") + ggtitle("Headline-specific tips effect") + ylim(-.7, .7) +
  coord_flip()

gg_slope <- g + theme_minimal() + theme(axis.text.y = element_blank())

ggboth <- plot_grid(gg_intercept, gg_slope, rel_widths = c(1.4, 1))
ggboth

# ggsave(ggboth, file = "headline_ranfx.pdf", width = 8, height = 7)

